
function [eqL,exitL,eqU,exitU] = EQ_bounds(gamma,B,X,xi,options,cmax,M)

% smallest equilibrium
diff = 1; exitL = 0;
eqL = zeros(5,1);
while diff >= 1e-2 && exitL <= 50
    x = BRall(gamma,B,X,eqL,xi,options,cmax,M)';
    diff = norm(x - eqL);
    eqL = x;
    exitL = exitL + 1;
end
    
if nargout > 2

    % largest equilibrium
    diff = 1; exitU = 0;
    eqU = repmat(cmax,5,1);
    while diff >= 1e-2 && exitU <= 50
        x = BRall(gamma,B,X,eqU,xi,options,cmax,M)';
        diff = norm(x - eqU);
        eqU = x;
        exitU = exitU + 1;
    end

end



